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ABSTRACT 

A model of the density distribution in the intergalactic medium, motivated 
by that found in numerical simulations, is used to demonstrate the effect of 
a clumpy IGM and discrete sources on the reionization of the universe. In 



an inhomogeneous universe reionization occurs outside-in, starting in voids 



and gradually penetrating into overdense regions. Reionization should not be 



sudden but gradual, with a continuous rise of the photon mean free path over 
a fair fraction of the Hubble time as the emissivity increases. We show that a 



hydrogen Gunn-Peterson trough should be present at z ~ 6 unless the emissivity 
increases with redshift at z > 4. However, the epoch of overlap of cosmological 



Q H - H ii regions could have occurred at a higher redshift if sources of low luminosity 

reionized the IGM; the Gunn-Peterson trough at z ~ 6 would then appear 
because even the most underdense voids have a large enough neutral fraction 
in ionization equilibrium to be optically thick to Lya photons. Cosmological 
H ii regions near the epoch of overlap can produce gaps of transmitted flux only 
^ ■ if luminous quasars contributed to the reionization, producing large H n regions. 

Despite the dumpiness of the matter distribution, recombinations do not 
increase the required emissivity of ionizing photons by a large factor during the 
reionization of hydrogen because the high density gas is not ionized until a late 
time. We show that the Hen reionization was most likely delayed relative to 
the hydrogen reionization, but was probably complete by z ~ 3 (the redshift 
where observations are available). The reported large optical depth fluctuations 
of He ii are not necessarily due to an incomplete He ii reionization, but can arise 
from a combination of IGM density fluctuations and variations in the intensity 
of the Hen ionizing background due to luminous QSO's. 



Subject headings: cosmology: theory — intergalactic medium — large-scale 
structure of universe — quasars: absorption lines 
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1. Introduction 

In the standard Big Bang model, the primordial gas becomes neutral at the 
recombination epoch, at redshift of z ~ 1100. The absence of the hydrogen Gunn- Peterson 
trough (which would be expected in high redshift objects if hydrogen were neutral) implies 
that most of the hydrogen in the intergalactic medium (hereafter, IGM) was highly ionized 
at redshifts z ^, 5 (Gunn & Peterson 1965; Bahcall & Salpeter 1965; Scheuer 1965; see 
Schneider, Schmidt, & Gunn 1991; Dey et al. 1998; Weymann et al. 1998; Fernandez- Soto, 
Lanzetta, & Yahil 1999, for the most recent evidence at the highest redshifts). The double 
ionization of helium may, however, take place later. Four quasars have been observed 
so far at the Hen Lya wavelength near a redshift z = 3 (Jakobsen et al. 1994, 1996; 
Davidsen, Kriss, & Zheng 1996; Hogan, Anderson, & Rugers 1997; Anderson et al. 1998; 
Heap et al. 1999), which show that a large fraction of the flux is absorbed (about 75% over 
2-5^-2^3); in addition, Heap et al. (1999) show that over 98% of the flux is absorbed in 
Q0302-003 at z > 3, except for some individual, well resolved gaps. Reimers et al. (1997) 
suggested that variations in the ratio of the He n to the H i optical depth provide evidence 
that the overlap of Hem regions was not yet complete at z ~ 3, but this is still controversial 
(Miralda-Escude 1998, Anderson et al. 1998; Heap et al. 1999). 

The IGM is well known to be highly inhomogeneous, owing to the non-linear collapse of 
structure. So far, most analytic treatments of the reionization process have only considered 
the inhomogeneity of the IGM in terms of a clumping factor that increases the effective 
recombination rate (e.g., Arons & Wingert 1972; Shapiro & Giroux 1987; Madau, Haardt, 
& Rees 1998). In this paper, we address the effects of inhomogeneity in more detail. We 
shall argue that another important factor is that reionization should occur gradually over 
the time interval during which the emissivity increases, and that denser gas will tend to 
be ionized at a later time than lower density gas. This can result in the ionization of 
most of the volume of the universe and the overlap of ionized regions, without the need 
to ionize very dense gas and therefore without the need to greatly increase the number of 
recombinations owing to the dumpiness of the IGM. 

The inhomogeneity of the gas distribution is also of primary importance to understand 
the absorption spectra of high redshift sources. The mean flux decrement of the photoionized 
intergalactic medium (or Lja forest) increases with redshift. Eventually, as the epoch of 
reionization is approached, the Lya forest should turn into a series of "clearings" , or "gaps" 
of transmitted flux appearing through the otherwise completely opaque Gunn-Peterson 
trough (Meiksin & Madau 1993; Miralda-Escude 1998), and the number of clearings should 
decline rapidly with increasing redshift. Quantitative predictions are needed to assess how 
this transition takes place in detail, depending on the spatial distribution of gas in the IGM 
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and on the luminosity function and lifetime of the sources of ionizing photons. 

In recent years, numerical simulations have led to considerable progress in the 
understanding of the density distribution and the ionization state of the IGM (Cen et al. 
1994; Hernquist et al. 1996; Miralda-Escude et al. 1996; Zhang, Anninos, & Norman 1995; 
Zhang et al. 1998). First attempts have also been made to incorporate ionizing sources 
and to investigate the corresponding evolution of the ionizing background during the epoch 
of reionization (Gnedin & Ostriker 1997). Here we use this knowledge of the density 
distribution of the IGM to make a simple model of the way reionization proceeds. 

In §2 we define the concept of the global recombination rate, and show how it determines 
the advance of reionization. In §3 the evolution of the intensity of the background is 
expressed in terms of the mean-free path, and we discuss the effects that determine the 
number of gaps of transmitted flux visible in Lya spectra as the epoch of reionization is 
approached. In §4 we study the conditions under which QSO's or star-forming galaxies can 
meet the requirements for reionization, and we discuss the question of the detectability of 
single H n regions around individual sources in the Lja spectra. §5 discusses a number of 
special aspects of helium reionization, and §6 presents our conclusions. Unless otherwise 
specified, we use the cosmological model H = 65kms _1 Mpc^ 1 , f2 = 0.3, A = 0.7, 
n h h 2 = 0.019. 

2. Reionization of an Inhomogeneous Intergalactic Medium 

The first discussion of the process of the reionization of hydrogen in the universe 
assumed an IGM with uniform density (Arons & Wingert 1972). It was immediately clear 
that the ionization would not take place homogeneously, but in individual cosmological 
H ii regions growing around the first sources of ionizing photons. The mean free path of 
ionizing photons in the neutral IGM, A, is very small [A = 0.28 (VL^h 2 /0.02)" 1 (1 + z)~ 3 Mpc 
at the threshold frequency, in proper units]. Unless the sources were very numerous and 
of low luminosity (with a number density as high as A~ 3 ), the thickness of the ionization 
fronts should be very small compared to the size of the Hn regions, and the IGM should 
be divided into H n regions around every source and neutral regions in between. As the 
mean emissivity in the universe grows, the Hn regions grow in size until they overlap. Two 
conditions need to be satisfied before the IGM is completely ionized: one photon for each 
atom in the IGM needs to be emitted, and the mean photon emission rate must exceed 
the recombination rate so that the hydrogen atoms are photoionized faster than they can 
recombine. If recombinations are never important, the Hn regions will never reach their 
Stromgren radius before overlapping. If recombinations are important, then the Hn regions 
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will fill an increasing fraction of the IGM as the emissivity rises and the recombination 
rate is reduced due to the expansion of the universe. In both of these two cases (and 
assuming that the size of the Hn regions is much smaller than the horizon), the model of 
the homogeneous IGM implies that there is a universal time when the last neutral regions 
are ionized. The mean free path of the ionizing photons and the intensity of the background 
increase very abruptly at this time. The reason for this sudden increase is that the universe 
is still very opaque to ionizing photons when only a small fraction of the atoms remain to 
be ionized, and then it becomes transparent in the very short time it takes to ionize these 
last atoms. 

Based on this model, the argument has been put forward that there should be a 
well-defined "redshift of reionization" at which this sudden rise of intensity takes place (e.g., 
Haiman & Loeb 1998 and references therein). This conclusion is, however, an artifact of 
considering a model of the IGM with uniform density. Furthermore, even in a homogeneous 
model, only the initial increase of the ionizing flux will occur on a timescale short compared 
to the Hubble time. The build-up of the final flux and the decrease of the mean opacity to 
values of order unity still takes a fair fraction of the Hubble time for any realistic redshift 
evolution of the emissivity of ionizing photons. 

We shall find in this paper that, when studying the problem of how the increasingly 
dense Lya forest turns into a Gunn-Peterson trough at high redshift, the assumption of an 
IGM with homogeneous density leads to a misleading picture of what actually happens. 
The inhomogeneous IGM must be considered right from the start, and it is then found that 
reionization must be a gradual process where the intensity of the background increases over 
a similar timescale as the emissivity of ionizing radiation. 

Let us therefore consider the volume-weighted probability distribution of the overdensity 
in the IGM, Py(A), where A = p/p, p is the gas density and p is the mean density of 
baryons. A crucial concept in our analysis will be that of the global recombination rate. 
Imagine that all the gas at density A < A ; is ionized, while the higher density gas is 
neutral. This assumption clearly cannot be exact, because the degree of ionization of 
gas of fixed density will depend on the local intensity of ionizing photons and the degree 
of self-shielding; nevertheless, we shall see that it provides a useful approximation to 
understand how reionization proceeds. Then, the mean number of recombinations per 
Hubble time that take place for each baryon in the universe is 



Here, R u is the ratio of the recombination rate for a homogeneous universe to the Hubble 
constant, given by R u (z) = 0.035 [Q b h(l - 3F/4)/0.025] (1 + zfl 2 /^ 2 for hydrogen 




(1) 
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(we use the high redshift approximation Q(z) ~ 1, and a recombination coefficient 
a = 4 x 10 -13 cm 3 sec -1 , for temperature T = 10 4 K; the factor 1 — 3Y/4 is the ratio 
of electron to baryon density, valid when helium is only once ionized). For Q = 0.3 
and h = 0.65, R u = 1 at z ~ 6. We neglect here the effect of the mean variation of 
the temperature with the density, which can affect the recombination rate due to the 
temperature dependence of the recombination coefficient, and the presence of collisional 
ionization. These effects become important only at low redshifts, when halos with a long 
cooling time for the hot gas are formed. 

Now, we consider how reionization should proceed in an inhomogeneous medium 
as the sources of ionizing photons appear. Most of the sources will be located in high 
density regions (galaxies), where stars or quasars can form. There is therefore a first phase 
of reionization where sources must ionize the dense gas in their host halos. This local 
absorption can simply be subtracted from the emissivity, so that the effective emissivity for 
reionizing the IGM counts only the emitted photons that are not absorbed in the halo where 
the source is located. Low luminosity sources may never be able to ionize the surrounding 
dense gas, and will then not contribute to the reionization. Thus, the period of reionization 
of the IGM starts when some sources have ionized their host halos and start emitting 
ionizing photons to the IGM (notice that the emission will then be anisotropic even for 
intrinsically isotropic sources, because the local absorption in the host halo will generally be 
anisotropic). The cosmological Hn regions will expand fastest along the directions of lowest 
gas density, both because fewer atoms need to be ionized per unit volume, and because 
fewer recombinations take place. As a simple example, if the gas density p g is constant 
along a narrow solid angle from the source (and for isotropic emission), then the distance 
d reached by the ionization front along different directions at a fixed time is d oc p g 1 ^ 3 
when recombinations are neglected, and d oc p g 2 / 3 when the recombination time is short 
compared to the age of the ionized region. 

The end of this period of reionization is reached when the ionized regions start 
overlapping, and a typical site in the IGM is illuminated by more than one source. Because 
the ionization fronts advance fastest along regions of low density, the overlap of H n regions 
will naturally occur first through the lowest density "tunnels" found between sources, and 
the denser gas with A > Ai will be left in neutral clumps. At the epoch of overlap the 
photon mean free path will be of the order of the separation between sources (i.e., the size of 
the Hii regions), and thereafter it will continue to increase as the characteristic density A ; 
increases. The dense gas clumps with A > A ; should then be the observed Lyman limit and 
damped absorption systems. The intensity of the background is limited by the absorption 
in these systems, until the photon mean free path increases up to the horizon scale. At this 
last stage, the universe becomes transparent and the intensity of the background is then 
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controlled by the redshifting of the radiation. 

Of course, the assumption that the stage of ionization of the gas is a function of the 
gas density only is a crude approximation. Some dense clumps will be ionized early if they 
are close to a luminous source, and the degree of self-shielding depends on the size of the 
clump in addition to its density. However, we shall see shortly that there is a simple reason 
why this approximation should still be a useful one. 

We will now write an equation governing the reionization of this inhomogeneous IGM. 
Let e(t) be the mean volume emissivity, measured in terms of the number of ionizing 
photons (with energy above 13.6 eV) emitted for each atom in the universe per Hubble 
time. We define also the ratio of the mean density of ionizing photons present in the cosmic 
background to the mean density of atoms, rij. Then we have 



where i^(Aj) is the fraction of mass with density A < A ; . The last term includes the losses 
of ionizing photons due to redshift. The constant (f> is defined as (p = [I v (vt)] [f^ dv [l v j u)} 1 , 
where I v is the background intensity per unit frequency and vt is the threshold frequency 
for ionization. Equation |2| reflects the statement that for every ionizing photon that is 
emitted, there must be either a photon added to the existing ionizing background, or a 
new atom being ionized for the first time, or an atom recombining to compensate for the 
absorption of an ionizing photon, or a photon being redshifted below the ionization edge of 
hydrogen (we neglect here the presence of helium) . 

In this paper we shall consider only the limiting case rij <C e, which is valid when the 
mean free path of the ionizing photons, Ai, is much smaller than the horizon, and provided 
also that A;/c is much smaller than the timescale over which e(t) increases. These conditions 
should generally be valid before the universe becomes transparent (unless the sources evolve 
synchronously over a timescale much shorter than the Hubble time). The terms containing 
nj in equation (0) can then be neglected, so 




In this approximation, the background intensity is given by 



and depends only on the instantaneous value of the emissivity (and not its history) because 
all the photons are absorbed shortly after being emitted. 




(2) 





(3) 
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Armed with equation (||), we can now see the reason why the basic approximation 
in our model (that gas with A < Aj is ionized, and gas with A > Aj is neutral) is 
reasonable. The crucial fact is that -R(Aj) is an increasing function of A i5 so that the global 
recombination rate is always dominated by the ionized gas with densities near Ajj we shall 
see in §2.2 that this is true for the expected density distribution in the IGM. Initially, when 
the first sources appear, the term dF M (Ai) / (H dt) dominates as gas is ionized for the 
first time. As previously mentioned, even in this regime the gas at low densities is ionized 
faster than the high density gas. But as the emissivity increases, R(A{) will become the 
dominant term at some point, as long as the photon mean free path remains smaller than 
the horizon (which we know does not occur until z ~ 2). At this point, the balance between 
the emissivity and the global recombination rate determines the density A ; up to which 
the gas is ionized as the emissivity rises. Most of the gas with A > A ; cannot possibly be 
ionized, because the recombination rate would then exceed the emissivity and the dense gas 
would rapidly become neutral. At the same time, most of the gas with A < A; needs to 
be ionized, because as we shall see in §3, the mean free path between regions with A > A ; 
increases rapidly as a function of Aj. Therefore, the typical photon emitted to the IGM will 
spend most of its trajectory in regions with A <C Aj before reaching a region with A ~ A ; 
and being absorbed, implying that if low-density regions were not ionized most photons 
would need to be absorbed there. In reality, the fraction of gas that is ionized will of course 
vary gradually over a range of densities around A i; and the effect of this gradual variation 
is what we neglect by approximating it as a sudden transition at density A;. 

This approximation is of course much less accurate before the overlap of H n regions, 
owing to the large fluctuations in the intensity of ionizing radiation. A better model is 
then obtained by incorporating the new variable Qi, the fraction of the volume in the IGM 
occupied by Hn regions. In the volume fraction Qi, some gas is still neutral because of its 
high density, whereas in the fraction 1 — Qi all the gas is neutral. Again, we can use the 
approximation that gas with overdensity A < Aj is ionized in the fraction Qi of the volume, 
and the higher gas density is neutral. Equation (|3]) is then modified to 

Notice that e/Q\ is simply the emissivity of sources averaged over the fraction Qi of the 
volume only. As before, there is a first epoch where the two terms involving the variation 
of Qi and Fm are dominant, and at some point recombinations start dominating globally 
to balance the emission. If Qi is still less than unity at this point, we can justify the use of 
equation (Jl]) for R(A{) in the following way. In every Hn region around individual sources, 
the radiation along every direction ionizes gas up to a point of high enough density where 
atoms recombine at the same rate at which the photons are arriving. If R(A{) is dominated 
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by the highest densities close to A i; most of the photons are absorbed in a dense region at 
the end of every beam. Once recombination dominates, these regions should have a typical 
density A ; which is still determined by a balance between the global rates of recombination 
and emission. There should of course be some gas with A > A, that is ionized in the 
clumps close to the luminous sources, but again the majority of photons will be absorbed 
at densities near A ; , where most of the gas clumps become self-shielding at a distance from 
the source equal to the mean radius of the Hn region. 

Equation ^ does not determine how Aj will increase with e before Qi reaches unity, 
because only one equation is provided for the two independent variables Qi and Ai once e(t) 
and F M (Ai) are known. The relative growth of Qi and Aj should depend on the luminosity 
function and spatial distribution of the sources. When the sources are very numerous, every 
void can be ionized by sources located at the edges of the void, and the overlap of ionized 
regions can occur through the thin walls separating voids, so Qi approaches unity when A; 
is the characteristic overdensity of these thin walls, Ai ~ 1. For rare and luminous sources, 
the value of Ai must be higher before Qi can be close to unity, so that the mean free path 
between the neutral clumps with A > Ai is equal to the mean separation between sources 
near the epoch of overlap. The early stages of reionization are characterized by the growth 
of Qi, while A; essentially determines the photon mean free path (which is of order the size 
of the Hn regions before the overlap). 

The end of reionization can be precisely defined as the moment when Qi = 1. In other 
words, after this moment there is no region in the IGM with low density (A <^ 1) which is 
not ionized. A word on what is meant by gradual reionization may be useful here. We note 
that Qi should reach unity rather fast at some universal redshift, because 1 — Qi will drop 
roughly exponentially with the mean number of sources seen from a random point, which 
is proportional to the cube of the mean free path. This implies an almost discontinuous 
change in the term in equation ||. However, the terms dF M (A^ / (H dt) + -R(A ; ) 
guarantee that A i; and therefore A; and rij, will not have a sudden increase at this epoch. 
Gradual reionization means that the logarithmic derivative dn j / (nj H dt) will not have 
a sharply peaked maximum over a very short time interval (t <C H" 1 ), but instead the 
increase of rij will take place over ~ a Hubble time. 

2.1. The Gas Density Distribution 

To examine more quantitatively the effects discussed above, we shall use a model of 
the gas density distribution in the IGM based on the results obtained from hydrodynamic 
numerical simulations, which were found to be in good agreement with observations of the 
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distribution of the transmitted flux in the Lya forest (Rauch et al. 1997). Here, we obtain 
a simple analytical fit to the numerical results of the density distribution in the IGM at 
different redshifts. 

The formula we use to fit the result of numerical simulations is motivated by a simple 
approximation of the evolution of the density in voids. In the limit of low densities, we 
assume the gas in voids to be expanding at a constant velocity, as should be the case if tidal 
forces are negligible. Then, the proper density in the void decreases as Aa -3 oc t~ 3 oc a~ 9 / 2 , 
so A a a~ 3 / 2 (where we assume £l(z) — 1). Assuming that A starts decreasing according 
to this law when the linear overdensity S reaches a fixed critical value, and using S oc a, we 
obtain A oc (— <5)~ 3 / 2 . The linear density S is assumed to be smoothed on the Jeans length 
of the photoionized gas. For Gaussian initial conditions, where Pv($) oc exp[— S 2 /(2a 2 )} 
(with a being the rms density fluctuation), this yields Py(A) oc exp (-CA- 4 / 3 ) A" 8 / 3 dA, 
where C is a constant. We shall use the following fitting formula, which approaches this 
limit for A 1, and is also applicable in the linear regime when |A — 1| <C 1, 

A' 13 dA . (6) 

When 5q <C 1 and Co = 1, the distribution approaches a Gaussian in A — 1 with dispersion 
8 . When S 3> 1, the distribution goes to the previous approximation for voids, with the 
peak at a density A ~ 5 3 ^ 2 . Thus, if the median density in voids evolves as A oc a~ 3 / 2 , 
then So should grow proportionally to a both in the linear and the highly non-linear regime. 
At A ^> 1 (and for So ^> 1), the distribution is a power-law, corresponding to power-law 
density profiles A oc r~ 3// ^ -1 ^ in collapsed objects (where r is the radius). 

We show in Figure 1 the gas density distribution of the L10 simulation described in 
Miralda-Escude et al. (1996), together with the fits obtained to equation fl6|), at the three 
redshifts z = 2, 3 and 4. The two functions shown at each redshift, for the simulations 
and the result of the fit, are APy(A) and A 2 Py(A), which are the volume-weighted 
and mass-weighted probability density of A per unit log A, respectively. This density 
distribution was previously shown in Figure 7 of Rauch et al. (1997) at z = 2 only, together 
with the result from a simulation by Hernquist et al. (1996). As seen in that paper, the 
density distribution in these two simulations is very similar, since they used models with 
very similar power spectra. The density distribution should depend mostly on a single 
parameter, the amplitude of the initial density fluctuations smoothed on the Jeans scale of 
the photoionized gas. Table 1 gives the values of the four parameters of the fits at each 
redshift. The parameter So is indeed approximately proportional to the scale factor. 

We have extrapolated the density distribution to higher redshifts by using 



iV(A) dA = A exp 



A- 2 / 3 - C c 
2 (25 /3) 2 
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z = 2 




log A 



Fig. 1. — Mass and volume- weighted differential density distribution of the IGM, at the 
three indicated redshits. The lines with noise are from the numerical simulation in MCOR, 
and the smooth lines are the fit we have obtained with the analytical model in equation (|6|). 
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Sq = 7.61/(1 + z), which reproduces the fits in Table 1 to better than 1%. The change of (3 
with redshift is more uncertain (the form of the high-density tail of the distribution in the 
simulations may also be affected by the resolution); we use (3 = 2.5 at z = 6, corresponding 
to an isothermal slope for high density halos. The parameters A and Co are then fixed by 
requiring the total volume and mass to be normalized to unity. 

We note here the possibility that a significant fraction of all baryons is contained 
in small halos with virial temperatures ~ 10 4 K, which collapsed to high densities before 
reionization, when the Jeans mass was very small (e.g., Abel & Mo 1998). These halos 
(which would not have been resolved in the simulations used for Figure 1) could survive for 
a long time after reionization occurs if neither star formation nor tidal and ram-pressure 
stripping due to mergers into larger structures are able to destroy them. In this case, 
the density distribution would be wider than in our model at high redshift, with a larger 
fraction of the baryons at very high densities. 



2.2. The Global Recombination Rate 

In Figures 2(a,b,c), the solid and dotted lines show the cumulative probability 
distribution of the gas density weighted by mass and volume, respectively, at redshifts 
z = 2,3,4 . These are obtained from equation (§), with the fit parameters in Table 1. 
Figure 2d shows the same cumulative distributions obtained at z — 6, according to the 
prescription described for choosing the parameters in equation (|6|); the values of these 
parameters at z = 6 are also given in Table 1. The dashed line shows the ratio R/R u as 
obtained from equation (fj). 

The model used in equation ([!]) sets the ionized fraction to unity for A < Aj, and 
to zero for A > A;. Clearly, the change in the mean ionization with A should in reality 
be more gradual. To test the sensitivity of the model to this assumption, we use also 
an alternative formula for the global recombination rate based on ionization equilibrium 
with a constant background intensity (notice that the change of the ionized fraction with 



Table 1 



Redshift 


A 


So 


P 


Co 


2 


0.406 


2.54 


2.23 


0.558 


3 


0.558 


1.89 


2.35 


0.599 



4 0.711 1.53 2.48 0.611 



6 



0.864 1.09 2.50 0.880 
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Fig. 2. — Solid and dotted lines are the cumulative density distribution weighted by mass and 
volume, respectively. The dashed line gives the global recombination rate R when all the gas 
is assumed to be neutral at overdensities greater than A, and ionized at lower overdensities. 
The dash-dot line is the global recombination rate when ionization equilibrium with a uniform 
ionizing background is assumed, with the neutral fraction being one half at overdensity A. 
The values of R u at the redshifts of the four panels, z — 2, 3, 4, 6, are R u = 0.32, 0.50, 0.71, 
1.18 for hydrogen, and R u = 1.89, 2.98, 4.19, 6.99 for helium (a gas temperature of 10 4 K 
has been assumed). 
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density should be steepened by self-shielding effects, and widened by fluctuations in the 
background intensity and a dispersion in the size of the clumps). The ionized fraction x is 
then given by R u x 2 A 2 = TA(1 — x), where T is the photoionization rate, and collisional 
ionization is neglected. Defining A; as the overdensity where x — 1/2 gives A; = 2T/R U , 
and x = A;/(4A) [(1 + 8A/A;) 1 / 2 — 1], and the global recombination rate is given by 
i?(A;) = R u / °° dA Py(A) x 2 A 2 . This is shown as the dash-dot line in Figure 2. The result 
is very similar to the case of the sudden change of the ionized fraction, except for A; ^ 1 
(where the model is not relevant anyway because even at the earliest stages of reionization, 
with Qi <C 1, A; should be of order unity or larger). 

As long as Py(A) is less steep than A -3 at large A, -R(Aj) should increase monotonically 
with Ai, as seen in Figure 2. This property of the global recombination rate is the main 
reason why high-density regions should generally be ionized at a later stage than low-density 
regions when there is a balance between global recombination and emissivity. If low-mass, 
dense halos formed before reionization and not resolved in the simulation we use were 
present, then the wider density distribution implied would cause a more rapidly increasing 
global recombination rate with Aj. 

The late ionization of the high-density gas implies that the dumpiness of the IGM does 
not necessarily increase the number of photons required to complete reionization. As an 
example, if overlap occurs at z = 6 and with A^ = 3, only about 70% of the baryons need to 
have been ionized, but 95% of the volume is ionized (Fig. 2d). The dashed curves also show 
that R/R n = 1 at A; ~ 3, so dumpiness in this case does not increase the net number of 
recombinations, and actually makes reionization easier by reducing the fraction of baryons 
that need to be ionized. In fact, if additional gas is in small-scale, high-density halos not 
resolved in the simulations, the required number of photons to complete reionization should 
be further decreased. Since R u ~ 1 at z = 6, only about one ionizing photon per baryon 
needs to be emitted to complete reionization. A low value of A ; at the epoch of overlap 
requires that the sources are numerous; reionization by more luminous sources would imply 
a higher A;, and therefore a greater number of recombinations. 



3. The Lyman Alpha Flux Decrement 

We now address the question of the flux decrement that should be observed to the 
blue of the Lja wavelength, as the redshift increases toward the epoch of reionization. The 
optical depth of a uniform, completely neutral IGM, tq, is extremely large, 



r = 2.6 x 10 5 [Q h h(l - y)/0.03] [H (l + z) 3/2 /H(z)] [(1 + z)/lf 2 



(7) 
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Because of this, a small neutral fraction is sufficient to yield a high enough optical depth to 
make the transmitted flux undetectable, so the Lya spectrum might already be completely 
blocked at redshifts lower than the epoch of overlap. 

3.1. The mean free path of the ionizing photons 

In the simple model where all the gas at densities A > Ai is still neutral, and all the 
gas at lower densities is ionized, a photon is absorbed whenever it enters a region with 
A > Aj, assuming that these regions are sufficiently large to be optically thick. The mean 
free path A; is then simply the mean length of regions with A < A; along random lines of 
sight. Of course, in reality there should be a gradual change of ionization, and differences 
in the geometrical shape of the structures will cause varying degrees of self-shielding at a 
given gas density, so photons will be absorbed over a range of densities. Nevertheless, since 
most of the recombinations take place at densities near A;, most of the photons will be 
absorbed by gas near this density. The spacing between successive contours at Ai should 
then provide a reasonable estimate of the mean free path. Notice that we assume that A; 
is larger than the mean free path for the case of a uniform IGM, A u = [n a (x(l — Fm)]" 1 , 
where n a is the density of atoms, a is a frequency-averaged photoionization cross section, 
and Fm is the fraction of mass in gas with A < Aj. If A u > Ai, then the neutral clumps are 
optically thin and the mean free path is close to A u independently of the density structure 
of the IGM. We also emphasize that Ai is the mean free path only for photons that escape 
into the IGM. Sources of photons will generally appear in regions of high density, and some 
fraction of the photons will be absorbed locally; these photons are not counted in the mean 
emissivity e in equations (0) and (|3|). 

Here, we shall use a very simple model for the mean free path of ionizing photons. We 
assume 

A i = A [l-F y (A i )]- 2/3 , (8) 

where Fy (Ai) is the fraction of the volume with A < Ai. In the limit of high densities, when 
1 — Fy <C 1, this is valid for any population of absorbers where the number density and 
shape of isolated density contours remains constant as Ai is varied: the fraction of volume 
filled by the high density regions is 1 — Fy, so their size is proportional to (1 — Fy) 1 / 3 , and 
the separation between them along a random line of sight is proportional to (1 — Fy)~ 2 / 3 . 
The shape of the absorbers in large-scale structure theories is actually highly complex, with 
structures varying from ellipsoidal to filamentary to sheet-like as Ai is reduced. However, 
we have found that the proportionality Aj oc (1 — Fy)~ 2 ^ 3 is still approximately obeyed in 
the results of numerical simulations with photoionized gas dynamics, for Aj ^ 1 (see e.g., 
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Fig. 3 of MCOR). 

In this paper, we shall use the scale X H = 60kms _1 , which reproduces the scales 
of the Lya forest structures in the simulation of the cold dark matter model with a 
cosmo logical constant at redshift z = 3 in MCOR (designated as L10 simulation in that 
paper). Although the scale Ao should vary with redshift and with the cosmological model, 
the quantity Ao H stays roughly constant in this simulation. In fact, XqH is basically 
determined by the Jeans length of photoionized gas at high redshift, when the amplitude of 
density fluctuations reaches non-linearity at the Jeans scale; it should then increase at lower 
redshifts as structures collapse on larger scales. This increase might not be fully reflected 
in numerical simulations like those in MCOR owing to the small size of the simulated box. 



3.2. Optical depth of the ionized IGM 

Our next step is to compute the optical depth to Lya scattering through the ionized 
IGM, 71, in a region of density A. This is equal to To A (eq. 0), times the neutral fraction, 
which we assume to be in ionization equilibrium. For a constant gas temperature, this 
optical depth is 

a A 2 

n = r - -=- , 9 

ocnj fj 

where a is the recombination coefficient (we assume full ionization and neglect the 
effect of double helium ionization on the electron density), a is the frequency-averaged 
photoionization cross section, and we have defined fj to be the ratio of the local photon 
density of the ionizing background at a given point in the IGM to the mean photon density 
rij, given by equation (f|). The local intensity fluctuates due to the discreteness of the 
sources that can illuminate a given point. Before the Hn regions around individual sources 
overlap, the intensity fluctuations are of course very large, but they are rapidly reduced 
after the mean free path increases to values larger than the typical separation between 
neighboring sources (Zuo 1992; Zuo & Phinney 1993). 

Substituting equations and @ into @, we obtain for the optical depth of an 
ionized region, 

-i = -o-^-(l-^) 2/3 ^. (10) 
crc if A e/j 

Notice that Fy depends on A;, determining the maximum densities that are ionized, while 
A is the density in the region yielding an optical depth T\. The mean emissivity e can 
be obtained from equation (|3|). The term dF M /(Hdt) in equation (|3|) obviously depends 
on the model for the emitting sources, but if the sources do not evolve on timescales 
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much shorter than the Hubble time, then &Fm/{H dt) ~ 1. Here we shall simply use the 
expression e = 1 + R; the important point is that e = R when recombination dominates, 
and e ~ 1 when most of the baryons are being ionized if recombinations are not important. 
Thus, we have 

^ Q h h(l-Y) H,{l + zfl 2 ( l + z \^ c (1-F y ) 2 / 3 A 2 _ A 2 
Ti 0.03 H(z) { 7 ) HX 1 + R fj'^fj' W 

We have used a value a = 2 x 10~ 18 cm~ 2 , which should be approximately valid for the 
spectra emitted by quasars or star-forming galaxies. This equation can also be expressed in 
terms of the recombination rate for a uniform IGM, R u , 

- T ° f 1 "^ 2 ' 3 1.14 — ^— 7— — — — r-y— — , (12) 



an e X (l + R)/R u ' HX (1 +R)/R U 

where n e is the electron density. 

The value of T u clS 8b function of A; is shown in Figure 3, at redshifts z = 2, 3 and 4. 
We have used cj (H\ Q ) = 5000, and the cosmological model mentioned in the introduction. 
As Ai increases, the larger mean free path of the photons results in an increasing intensity 
of the ionizing background, and that decreases the optical depth of the ionized regions. 
The variation of r n (A{) with redshift can be understood as follows: ignoring the variation 
of the density distribution Fy with redshift, and for large A i; we have R ^> 1, and r u is 
independent of redshift from equation ([12]). For small A;, r u increases with redshift since 
i?< 1. 

If the IGM were homogeneous, the mean flux decrement due to Lya scattering would 
simply be e~ Tu , but for the real IGM the dependence of the flux decrement on r u is of course 
more complicated. The mean flux decrement of the Lya forest obtained from the simulation 
in MCOR is shown as a function of r u in Figure 4, at redshifts z = 2,3,4. These curves 
were computed directly from the simulation, and therefore include the effects of thermal 
broadening and peculiar velocity; they were shown also in Fig. 19 in MCOR (although with 
a different scaling in the horizontal axis). Each curve starts at the point where the flux 
decrement is equal to the observed value (see Press, Rybicki, & Schneider 1993; Rauch et 
al. 1997). Therefore, in the MCOR model, if all the baryons were spread out uniformly in 
the IGM keeping the intensity of the ionizing background fixed, the optical depth to Lya 
scattering would be r u = (0.3, 1, 4) at z — (2, 3, 4). The curves show how the flux decrement 
would increase at a fixed redshift if the intensity of the background was reduced, increasing 
the value of r u . At fixed r u , the flux decrement increases slightly with redshift. The reason 
is that voids become more underdense with time, and therefore the amount of flux that is 
transmitted through the most underdense regions for high r u decreases. 
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Fig. 3. — Optical depth r u that a uniform IGM would yield for the intensity of the ionizing 
background derived when the emissivity and the mean free path are given by the global 
recombination rate and our mean free path model, as a function of the overdensity to which 
the inhomogeneous IGM is ionized. 




Fig. 4. — Mean flux decrement as a function of the optical depth that a uniform IGM with 
the same mean density and intensity of the ionizing background would have. 
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3.3. What is the redshift where the Gunn-Peterson trough is reached? 

As sources at progressively higher redshift are being discovered, one of the most 
interesting questions that arise is: what is the redshift at which a complete Gunn-Peterson 
trough will first appear? It is often assumed that the Gunn-Peterson trough should appear 
when the epoch of reionization is reached. However, a small neutral fraction suffices to make 
the IGM opaque to Lja photons everywhere, given the enormous value of the optical depth 
for a neutral medium (eq. [?]). In ionization equilibrium with a uniform background, the 
most underdense voids will have the lowest optical depths. For the gas density distribution 
used in §2, the lowest densities are A ~ 0.1 at z = 2 and A ~ 0.25 at z — 6. Since the 
optical depth is proportional to A 2 , the transmitted flux ought to decrease very rapidly to 
negligible levels when r u reaches the inverse square of the lowest underdensities in voids, as 
the curves in Figure 4 show (the peculiar velocity effect in voids helps to allow for some 
transmitted flux at slightly higher values of r u ). 

Figure 5 shows the same flux decrement as in Figure 4, plotted as a function of A ; . 
The large filled squares are the values of the observed transmitted flux, according to Rauch 
et al. (1997). The figure therefore provides a prediction for A; as a function of redshift, 
given our assumed model for the density distribution: the gas in the Lyman limit systems 
at column densities ~ 10 18 cm -2 , where the transition from ionized to neutral gas takes 
place due to self-shielding, should have typical overdensities of order Ai ~ (700, 300, 100) 
at z — (2,3,4). Notice, however, that the approximation of neglecting redshift effects in 
equation (fj) is starting to fail at z = 2, so A ; should be somewhat underestimated. 

It is clear from Figure 5 that, if the trend of decreasing A ; with redshift continues at 
z > 4, then the Gunn-Peterson trough should appear at z ~ 6. There are two main reasons 
for the increase of the flux decrement with redshift. The most important is the decrease of 
the density Aj to which the universe is ionized. The second one is that the flux decrement 
increases for fixed A ; owing to the lower underdensities of voids at high redshift; thus, at 
z ~ 6 the transmitted flux should drop to 1% at A ; ~ 15. 

Figures 2(a,b,c) show that the global recombination rate corresponding to the densities 
A; of Lyman limit systems mentioned above, at z — (2,3,4), are R/R u = (40,15,6), 
corresponding to an emissivity of e = R = (14, 8.1, 4.3) photons per baryon per Hubble time. 
If the trend of decreasing emissivity with redshift continues, clearly the Gunn-Peterson 
trough will be reached at z ~ 6. In fact, in order to have more than a few percent of 
transmitted flux at z — 6, according to Figures 5 and 2d, we need A; ^ 20 and e > 2.5. 
Since we defined e as the number of photons emitted per Hubble time, then as long as the 
emissivity per physical time is not larger at z = 6 compared to z — 4, the transmitted flux 
at z = 6 should be less than ~ 3%. We emphasize that the reason is not because the epoch 
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Fig. 5. — Mean flux decrement as a function of the overdensity up to which the gas is ionized. 
Solid squares indicate the observed value of the mean flux decrement according to Rauch et 
al. (1997) at z = 2,3,4. 
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of reionization is reached at this redshift, but because even the most underdense voids have 
a large enough density of neutral hydrogen to scatter essentially all the Lja photons from 
a background source. 

Nevertheless, the redshift at which the Gunn-Peterson trough is first seen might be 
higher if more sources were present at high redshift. For example, if the emissivity remained 
at e ~ 5, then A ; ~ 100 at z = 6 and the flux decrement would not decline very fast 
at z > 4. In fact, because less mass has collapsed to high densities at high redshift, the 
value of A; (and therefore the flux decrement) becomes extremely sensitive to R (Fig. 2d). 
However, our model of the density distribution does not take into account that dense gas 
could be present in low-mass halos that collapsed before the IGM was ionized, on scales 
below the Jeans scale of the photoionized gas. 

The observations of the highest redshift sources (Dey et al. 1998, Weymann et al. 1998) 
suggest so far that the flux decrement rapidly declines, indicating that the decrease of A; 
with redshift probably continues along the same trend shown in Figure 5 at higher redshifts. 

3.4. The mean free path and the abundance of Lyman limit systems 

Figure 6 shows the mean free path Aj as a function of Aj that we derive in our model. 
For the values of Ai obtained previously to match the observed flux decrement, the mean 
free path at z = (2,3,4) is H A, = (50,33,18) x 10 3 kms _1 , which is in good agreement 
with the number of observed Lyman limit systems per unit redshift, ~ 2[(1 + z)/4] L5 (e.g., 
Storrie-Lombardi et al. 1996). We notice that in the MCOR simulation, the number of 
Lyman limit systems was much less than observed. The reason we obtain a good match here 
is because we have fitted the high-density tail of the density distribution with a power-law 
that is close to an isothermal slope for gas distributed in halos, and we have also used a 
power-law dependence of Aj on A ; that assumes the presence of halos with density cusps. 
The numerical simulation in MCOR contains instead halos of gas with core sizes limited by 
the resolution, which explains why the number of Lyman limit systems in the simulations 
was smaller than in our model here. 



3.5. Gaps in the Gunn-Peterson trough due to individual Hn regions 

Is it possible that some of the last gaps of transmitted flux appearing in Lja spectra 
will correspond to individual H n regions around the sources that reionized the IGM before 
the epoch of overlap? As pointed out in Miralda-Escude (1998, hereafter M98), there are 
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Fig. 6. — Mean free path of ionizing photons as a function of the overdensity up to which 
the gas is ionized, according to equation @, at the indicated redshifts. 
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two conditions that need to be satisfied for an ionized region to transmit flux. The first 
condition is the one we have also discussed here: the optical depth through the ionized 
region should not be much higher than unity. Given the discussion at the end of §2, the 
formalism developed here to compute the typical value of n from equation ([TI|) should be 
equally valid before the epoch of overlap, in Hn regions around individual sources, taking 
into account the following differences: (a) The mean emissivity is now substituted by e/Qi, 
averaging only over the ionized volume of the universe, (b) The mean free path Aj should 
be identified with the radius of the H n region, (c) The fluctuations of the intensity should 
be much higher than after the epoch of overlap, so when the line of sight passes near a 
luminous source 73 may be decreased owing to a high value of fj in equation flll]) . 

From Figure 4, we see that the transmitted flux rapidly becomes very small when 
the optical depth of a uniform IGM with the same intensity of the ionizing background 
is r u ^ 50. For our model of the mean free path, this corresponds to Aj <^ 20, and 
H(z)\i ^ 4000 km s -1 near z = 6. This condition is equally applicable to individual 
H 11 regions: in order for an H 11 region to be visible in the Lya spectrum as a gap in 
the Gunn-Peterson trough, the mean intensity of radiation in the Hn region should be 
large enough to make r u <^ 50, implying (for z = 6) that the gas must be ionized up to 
densities Ai ~ 20 in order for the ionized region to be at least as big as the mean free path 
A; ~ (4000kms~ 1 )/iJ(z) given by Figure 6. The spectrum of the line of sight through an 
Hn region of at least this size will then show a gap in the transmitted flux when voids are 
crossed with the average ionizing intensity given by the luminosity of the source and the 
size of the ionized region. Smaller H 11 regions can be seen only rarely if the line of sight 
passes through a void close to the source (so that fj 1). 

The second condition for an isolated Hn region to transmit any Lya photons arises 
from the fact that if a large fraction of the IGM is still neutral, then the damped absorption 
profiles of the neutral gas in front and behind the Hn region will overlap, scattering the 
remaining photons that could have crossed an underdense, ionized region. It was shown in 
M98 that an Hn region must have a proper size larger than lh^l^l^hil — K)/0.03] Mpc to 
prevent this overlap of damping wings, and that the optical depth in an H 11 region of this 
size is r u = 475(T/10 4 K)" a7 e" 1 (see eq. 4 in M98; the quantity N r in M98 is 2n h /3n e times 
our e, where n e /n^ = 1 — Y/2), independently of the redshift or any other parameters. In 
general, this condition is not as restrictive as the first one, r u <^ 50, needed for the ionized 
region itself to have a low enough optical depth, if we assume e~R~2&tz = G. However, 
for small H 11 regions this second condition implies that even if the gap of transmitted flux 
could be visible due to fj 1, the flux might still be removed due to the damping wings. 

We therefore conclude that whether Lya photons can be transmitted before the overlap 
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of H ii regions is complete depends mainly on the typical luminosity of the sources of 
ionizing photons. We shall discuss plausible sizes of the H n regions in §4. 

This section has presented a general formalism to calculate the evolution of r u 
with redshift, the epoch when the Lya forest should become a complete Gunn-Peterson 
trough, and the condition on the size of the individual Hn regions required for Lja flux 
to be transmitted through them. However, our numerical results depend on the density 
distribution Py(A) and the scale of the density structures in the IGM, Ao, which we 
have taken as constant in reshift space here. These quantities will need to be calculated 
separately for every cosmo logical model as a function of redshift. 

4. Hydrogen Reionization 

4.1. Sources of Ionizing Radiation : QSO's and Galaxies 

Among known, observed objects, QSO's and galaxies are the two candidates for the 
sources of ionizing photons that reionized the IGM. At z ~ 3 the total number of UV 
photons emitted by galaxies longward of the Lyman break outnumbers those produced by 
optically bright QSO's by a factor of approximately ten to thirty (Boyle & Terlevich 1998; 
Haehnelt, Natarajan, & Rees 1998). For photons capable of ionizing hydrogen, the balance 
is less certain. Observed bright, high-redshift QSO's show no Lyman break, whereas 
galaxies have a strong intrinsic Lyman break arising from the stellar atmospheres (of a 
factor ~ 5 for the usual mass function of stars observed in galaxies), and a very uncertain 
escape fraction of photons beyond the Lyman edge. For starbursts at low redshift, the 
escape fraction is reported to be of order 10% (Leitherer et al. 1995; Hurwitz, Jelinski & Van 
Dyke Dixon 1997). A moderately higher value would therefore be required for UV emission 
by star-forming galaxies to rival that of QSO's at z = 3 in terms of ionizing photons. 

The solid curves in Figure 7 show a fit to the luminosity function of QSO's at z = 4 
obtained by Pei (1995): d$ = $*/((L/L*)^ + (L/L*)h) dL/L*, with = 8.9 x 1CT 7 Mpc~ 3 , 

= 7.8 x 10 45 ergs _1 (B-band) and (3 h = 3.52. The luminosity function is plotted as 
follows: on the horizontal axis, instead of the luminosity L, we use the variable 
-Rhii = [3Aphot/(47mH)] 1,/3 , where iV p h t is the total number of ionizing photons emitted by 
the source for an assumed lifetime £q, tir is the comoving number density of hydrogen, and 
i?Hii is the comoving radius of the H n region that a source of luminosity L would produce 
in a completely neutral homogeneous medium if there were no recombinations. We obtain 
the number of photons from the relation A p hot = LtQ/E[ on , where E\ OIi is the mean energy 
per photon. We have assumed that the QSO's emit twice the luminosity in ionizing photons 



-25 - 



as in the B-band (a reasonable approximation for the typical observed quasar spectra), 
and have used E ion = 20 eV. Also shown in the upper horizontal axis is the corresponding 
expansion velocity at z — 4, H(z = 4) _R H n/(l + 4). 

On the vertical axis, the mean comoving volume emissivity of QSOs per unit natural 
logarithm of Run is plotted, in terms of the number of photons emitted per hydrogen atom 
per Hubble time. The curves in Figure 7 are an "energy function" (or the distribution 
of total energy emitted by sources), rather than a luminosity function. The energy 
function depends of course on the assumed lifetime, and we have plotted it for two values: 
£q = 3 x 10 7 years and £q = 8 x 10 8 years. For other values of £q, the curves simply shift 
horizontally by a factor proportional to tq 3 . The time 8 x 10 8 years is equal to half the age 
of the universe at z = 4; this can be considered as an upper limit for tq, because the number 
of quasars of given luminosity should be rapidly increasing with time at high redshift. All 
values are computed for our adopted cosmological model mentioned in the introduction. 

The faint end slope of the luminosity function is still rather uncertain, so we have 
plotted three curves for the values (3\ = (1.5, 1.65, 1.8), where f3\ = 1.64 is the value preferred 
by Pei. The value of the lowest luminosity for which the abundance of quasars has been 
determined observationally at z — 4, 10 46 ergs _1 , is indicated with an arrow. For lower 
luminosities, the model assumes that the luminosity function has the same shape as at 
lower redshift. 

The diamond symbols in Figure 7 show in the same way the mean comoving volume 
emissivity of UV-emitting galaxies. We use a Schechter function fit to the luminosity 
function at z = 3 (Dickinson et al. 1998), d$ = (L/L*)^ exp (— L/L*) dL/L*, with 

= 1.6 x l(T 3 Mpc~ 3 , L, = 1.4 x lO^ergs" 1 (defined as vL w at 1500 A), and (3 = -1.5. 
The luminosity function of UV-emitting galaxies does not evolve strongly between z = 3 
and z = 4 (Steidel et al. 1998). The ionizing luminosity is assumed to be a fraction 0.06 of 
vL v at 1500A, due to an intrinsic Lyman break of a factor 5, and assuming a fraction of 30 
% of emitted ionizing photons to escape into the IGM. The age of the galaxies is set equal 
to 8 x 10 8 years, or half the age of the universe at z = 4. 

The integrals under the curves of the energy functions give us the total emissivity, in 
terms of number of emitted photons per hydrogen atom per Hubble time. For both QSO's 
and galaxies, their total emissivity at z = 4 is e ~ 5. The radii of the Hn regions when the 
lifetime is set to t s = (eiJ) -1 have a special significance: over the time (eH)' 1 , one ionizing 
photon is emitted for each atom in the universe, and therefore the Hn regions around each 
source are just large enough for overlapping. Therefore, for this lifetime, the variable in the 
horizontal axis is equal to the radius that the ionized regions around a source of a given 
luminosity must have grown to at the epoch of overlap, when reionization ends. The curve 
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Fig. 7.— The solid curves show the mean comoving volume emissivity from QSOs in 
terms of the number of ionizing photons emitted per hydrogen atom per Hubble time as a 
function of the comoving radius of the reionized region created by a quasar over a lifetime 
£q, if recombinations are ignored. The QSO luminosity function at z = 4 of Pei (1995) 
is used. The three curves in each set are for faint end slopes f3\ = (1.5,1.65,1.8). The 
open diamond symbols show the mean comoving volume emissivity of star-forming galaxies 
at z = 4 calculated from the luminosity fuction according to Dickinson et al. (1998). An 
intrinsic Lyman break of a factor five and an escape fraction of ionizing photons of 30 percent 
is used to extrapolate from the flux at 1500 A and a lifetime £q = 8 x 10 8 yr is used. The 
arrows denote the radii which correspond to a 10 46 ergs _1 QSO. The dashed lines show the 
global number of recombinations taking place per Hubble time and per hydrogen atom as 
a function of the mean free path (equated to the H n region radius in the bottom axis) at 
redshifts z = 4 and z — 6. The recombination rate was calculated as described in §3 for the 
density distribution shown in Fig. 1. The dotted curves give the mass fraction required to 
ionize up to the overdensity corresponding to that mean free path, times a factor 3. 
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corresponding to this lifetime would be shifted to the left from the curve for ig = 8 x 10 s 
yr by a factor (5/3) 1 / 3 = 1.19. 

The requirement for having completed the reionization with a given population of 
sources can now be easily visualized from Figure 7. We shall consider first the case of 
long-lived sources that remain active with an approximately constant luminosity over the 
epoch of reionization; the differences in the case of short-lived sources will be discussed 
later. In the previous sections we computed the global recombination rate in terms of the 
maximum overdensity to which the universe is reionized, A;, which in turn can be expressed 
in terms of the mean free path A; (using equations [ED], ||, and M). The dashed lines in 
Figure 7 give the number of recombinations taking place per Hubble time as a function of 
Ai, equated here to the radius of the Hn region in the bottom axis. The dotted curves give 
the mass fraction Fm required to ionize the IGM up to the overdensity corresponding to 
the mean free path Aj, times a factor 3. These curves are shown at z = 4 (thick lines) and 
z = 6 (thin lines). At high redshift, the emissivity from quasars could of course be different, 
depending on how they evolve. At the epoch of overlap, the radii of the Hn regions must 
have grown to the values given in Figure 7 for the lifetime t$ = (eH)~ l . The first condition 
necessary to reach the epoch of overlap is that et^H > F M , where t R is the time during 
which sources have been producing the emissivity e, and etnH is the total number of 
photons emitted. Choosing tR = (3//) _1 (i.e., half the age of the universe at the epoch of 
overlap), this condition becomes e > 3F M , roughly equivalent to requiring that the solid line 
for a lifetime ts should be above the dotted line in Figure 7. The second condition is that 
the global recombination rate when the mean free path is equal to the size of the H n regions 
that must be reached for overlap is not greater than the emissivity. This condition roughly 
implies that the solid line must be above the dashed line. 

As the emissivity of the sources increases, the epoch of overlap is reached when the 
mean emissivity of the sources rises above both the dashed and dotted lines. This shows 
that the importance of recombinations, and the required emissivity to reionize the IGM, 
increases not only with redshift but also with the luminosity of the sources. 

For sources with luminosities typical of high-redshift starburst galaxies, the fact 
that the dashed line is below the dotted line implies that recombinations should not be 
very important in the case of late reionization by galaxies. Even if very luminous QSO's 
reionized the universe, recombinations increase the required number of emitted photons by 
a moderate factor only, because the IGM does not need to be ionized up to a very high 
overdensity to increase the mean free path up to the size of the H n regions produced by 
these sources. Recombinations are of course more important as the redshift of reionization 
is increased. 
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Now, we consider the case of short-lived sources, with lifetimes much less than tf>. 
As the reionization proceeds, the IGM should contain active H n regions, which continue 
to expand as more radiation is emitted by the source, and fossil Hn regions, which start 
recombining after their central source turns off, but are still mostly ionized. This implies 
that in the case of short-lived sources there are two different epochs of overlap, a first one 
for the overlap of fossil H n regions, and a second one for the overlap of active H n regions. 
The first epoch of overlap of fossil Hn regions is determined by the energy function of 
sources. This occurs when the emissivity rises above the dashed and dotted lines at the 
value of the mean free path equal to the radius of fossil Hn regions, determined by the 
lifetime £q. The second epoch of overlap is determined by the luminosity function, and it 
generally occurs significantly later and requires a higher emissivity, determined by the same 
condition as for long-lived sources obtained when the solid lines are shifted to the right 
to the lifetime t$ = (eH)' 1 . Between these two epochs of overlap, the emitted ionizing 
photons are invested not just in ionizing again the baryons that have recombined in fossil 
Hn regions, but in ionizing additional gas in more overdense regions in order to increase 
the mean free path to the separation between active H n regions. The shorter the lifetime 
of the sources, the more widely separated the two epochs of overlap become. 

In the example of a lifetime tQ = 3 x 10 7 years, fossil Hn regions overlap when a typical 
quasar with luminosity 10 46 ergs _1 has produced an Hn region of comoving radius 25Mpc. 
The overlap of active H n regions would be reached only when the radius of H n regions 
around a quasar of this luminosity increases up to a radius of 70 Mpc, requiring a higher 
emissivity to balance the global recombination rate. 

It is clear from Figure 7 that for both QSO's and galaxies, the emissivity is sufficient 
to have reionized the universe by z = 4. If we use the approximation e = R, our model 
predicts that at z = 4, A ; = 100 and R ~ QR U ~ 4 (see Figs. 2c, 5), so the combined 
emissivity of quasars and galaxies can actually be slightly higher than what is needed to 
match the observed flux decrement. At higher redshifts, the dashed line indicating the 
recombination rate increases and also becomes shallower, and the epoch of reionization is 
reached when the dashed line is equal to the emissivity at the mean separation between the 
sources. If the emissivity from low-luminosity galaxies remains at a constant level, this may 
not happen until z ~ 10 or higher. To make further progress, a determination of the escape 
fraction of ionizing photons from high-redshift galaxies and of the number densities of faint 
quasars and galaxies at high redshift will be required. 
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4.2. The Size of Hn Regions at the Epoch of Reionization 

We now come back to the question of the possible presence of gaps of transmitted flux 
due to individual H n regions in Lja spectra of high-redshift sources before the hydrogen 
reionization is completed. In §3 we showed that these gaps should generally be present if 
t u 50, which for our value of the mean free path corresponds to a size of the H n region of 
if Ai ^ 4000 km s" 1 , or (1 + z)X\ ^ 40 Mpc at z ~ 6; smaller Hn regions can also produce 
gaps when fj > 1. At higher redshifts, the required size of the Hn regions to produce a 
gap of transmitted flux increases, both due to the lower underdensities of voids, and to the 
longer mean free path at fixed r u (see Figures 3, 5 and 6). 

Figure 7 shows that the typical QSO's that dominate the emissivity at z = 4 are able 
to create Hn regions large enough to produce observable gaps in Lya spectra. In the case 
of short-lived QSO's, the Hn regions would be smaller at the first epoch of overlap of fossil 
regions. However, when the active Hn regions overlap, their size obviously depends only on 
the observed luminosity function of active sources, and not on their lifetimes. Star forming 
galaxies do not reach the required size, and therefore their associated Hn regions can only 
be visible in cases where fj is large. We notice here that the conclusion in M98 that these 
gaps of transmitted flux would probably not be present at redshifts above the epoch of 
overlap was due to assuming that the high-redshift sources responsible for reionization 
would be of low- luminosity compared to the known QSO's at lower redshift. 

The size of H n regions at overlap has strong implications for the effect of reionization 
on the fluctuations of the CMB background (e.g. Gruzinov & Hu 1998). Unfortunately the 
constraints are extremely weak. The typical size could range from less than 1 Mpc if faint 
star-forming galaxies were the ionizing sources to several tens of Mpc if reionization were 
caused by bright QSO's. 



5. Helium Reionization 

The double reionization of helium can be treated following the same method we have 
used for hydrogen. The Hem recombination rate is 5.5 times faster than that of hydrogen 
(this is increased to 5.9 times faster when including the increase in the electron density 
due to the ionization of helium), i? H ein,u(^) = 0.21 [Q h h{l - Y/2)/0.027] (1 + zf/ 2 /Q}J 2 . 
If the number of photons emitted above the Hen ionization threshold is not higher than 
5.9nHe/^H — 0.46 times the number of photons emitted above the hydrogen ionzation edge 
(which is the case for all existing candidates of the sources of the ionizing background), 
and if recombinations are the dominant balance to the emissivity, then the Hen should be 
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reionized at a later epoch than the hydrogen. 

Equations ([TJ), (|3|), (H), and (W) are equally applicable for Hen , once we define the 
new quantities en e and nj t u e as the emissivity per Hubble time and the number density of 
photons above the He n ionization edge for each helium atom, and An e as the density below 
which helium starts being doubly ionized. Equation (|i~2"D is also equally valid, except that 
the numerical factor changes from 1.14 to 1.23 (due to the ratio of the number of hydrogen 
atoms to the total number of hydrogen and helium atoms that is included in this equation 
for the case of hydrogen); the basic reason is that the ratio of the integrated optical depths 
to Lyo; line scattering and to continuum ionizing radiation is always the same for a Hn or a 
He in region. 

Figure 8 shows r u versus A Hc , and Figure 9 shows the flux decrement, at redshifts 
z = 2, 3, 4 (the latter is obtained directly from the numerical simulation in MCOR [see their 
Fig. 19], and it includes the effect of thermal broadening of helium). The symbols in Figure 
9 show the observed mean flux decrement in the following objects: triangle, Q0302 — 003 at 
z = 3.15 (Heap et al. 1999); square, Q0302 - 003 at z = 2.82 (Heap et al. 1999); hexagon, 
HS1700+64 at z = 2.4 (Davidsen et al. 1996). The value of An e indicated by each point in 
Figure 9 gives the overdensity of the gas up to which helium should be doubly ionized, in 
order to reproduce the flux decrement in our model. The implied mean free path of the 
Hen -ionizing photons can then be read off Figure 6. 

The results of the observations of the Hen Lya spectra can be summarized as follows: 
a flux decrement has been clearly detected in four quasars, which increases rapidly with 
redshift over the range 2.2 < z < 3.2. At the upper end of this range, a very small fraction 
of the flux is transmitted over most of the spectrum, but there are gaps (with widths 
~ 1000 km s -1 ) where the fraction of transmitted flux is high (Reimers et al. 1997, Anderson 
et al. 1998; Heap et al. 1999). Some of the wide gaps can be associated with the proximity 
effect of the observed sources; however, other gaps are far from the sources. According to 
Heap et al. (1999), outside of these gaps a very small fraction of transmitted flux of ~ 1% 
is detected up to the highest redshift observed (corresponding to the triangle in Fig. 9). 

Figure 10 is analogous to Figure 7, showing the emissivity of He n -ionizing photons 
per helium atom from quasars as a function of the comoving radii of their Hem regions, 
at z — 3, plotted for two assumed lifetimes: 3 x 10 7 and 10 9 years (the latter is ~ half 
the age of the universe at z — 3). The top horizontal axis gives the Hubble velocity 
corresponding to the radius of the Hem region at z — 3. We use again the parameterization 
of the QSO luminosity function given by Pei (1995). As before the dashed line shows the 
global recombination rate per helium atom. We have assumed a spectral index of 7 = 1.8 
(f v oc i/~ 7 ) to extrapolate from the hydrogen to the helium Lyman continuum (Zheng et al. 




8. — Same as Fig. 3, for the case of helium. 
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Fig. 9. — Mean flux decrement as a function of the overdensity up to which the gas is ionized, 
for the case of helium. Symbols indicate the observed mean flux decrements in the objects 
specified in the text. 
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1997). 

According to our model, the total emissivity of QSO's is sufficient to completely 
reionize Hen before z — 3, since the emissivity is above the global recombination rate for 
the lifetime ts = (euctl)' 1 ~ 1.7 x 10 8 yr. The mean free path should be given by the radius 
where the emissivity equals the recombination rate, A# e ~ 7000kms~ 1 , corresponding to 
An e ~ 50 according to Figure 6. This then implies a flux decrement ~ 0.8 from Figure 9; in 
addition, since the mean free path should already be larger than the separation between the 
quasars dominating the emissivity, fluctuations in the intensity of the background should be 
relatively small. These predictions seem to fit the observational results at a slightly lower 
redshift, z ~ 2.5. On the other hand, the observations of a much larger flux decrement with 
gaps of transmitted flux at z ~ 3 suggest a lower emissivity than in the model by Pei used 
in Figure 10. 

A model for the sources of He n -ionizing radiation that reproduces the observations 
at z ~ 3 should have a population of low-luminosity sources that has ionized the He n up 
to overdensities A He — 12, to reproduce the transmitted flux reported by Heap et al. 
(1999) outside the large gaps. From Figures 6 and 10, this implies a mean free path 
Arc — 1500 km s -1 and an emissivity eHe ~ 4 from these low-luminosity sources. Luminous 
quasars should also have a similar emissivity ene, placing it slightly below the dashed line in 
Figure 10, allowing for the presence of numerous gaps of transmitted flux caused by regions 
that are more intensely ionized by luminous quasars. The presence of sources with a wide 
range of luminosities can allow for large fluctuations in the intensity of the background due 
to the most luminous sources, well after reionization has been completed by less luminous 
sources. 

We notice here that, in the case of helium, recombinations are much more important 
than the need to ionize Hen for the first time, for luminous sources (i.e., the dashed line 
is well above the dotted line in Fig. 10). Because the recombination rate grows rather 
steeply with the mean free path, low-luminosity sources are very likely to have ionized 
the He ii entirely at z > 3. If the emissivity from sources with the abundance of galaxies 
in Figure 7 is more than 10% of the emissivity from luminous quasars, then galaxies 
should have reionized the Hen first, since the recombination rates at the mean separation 
between galaxies and luminous quasars differ by a factor ~ 10. It is quite plausible that 
the emissivity from galaxies reaches at least this level, because the cooling radiation from 
supernova remnants and hot gas in the halo produces photons of the right energy for the 
ionization of helium. Therefore, unless luminous quasars have always dominated by a large 
factor the emissivity of He n -ionizing photons, the He n should have been ionized by lower 
luminosity sources and the epoch of overlap of the He in regions should be earlier than 
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Fig. 10. — Same as Figure 7, for the case of helium, at z — 3. The emissivity from quasars 
is shown here in terms of number of photons emitted per Hubble time and per helium atom, 
as a function of the radius of the He in region that a quasar of lifetime £q would create in 
the absence of any recombinations. 
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z = 3. 

Thus, our interpretation of the gaps of transmitted flux in the He n Lja spectra is 
that they are due to an ionization effect of a luminous source, rather than simply a low 
density region in the IGM (which is supported by various tests presented by Heap et al. 
1999), but that they are not truly Hem regions before the epoch of overlap. Instead, 
the luminous quasars are only increasing the background intensity in a large region 
around them, and ionizing the dense gas in the absorption systems to greater depth. To 
illustrate this point, we can calculate the degree of ionization of the Hen when An e = 12, 
and therefore the transmitted flux in the Hen Lya spectrum is only 1% according to 
Figure 9. As mentioned above, this requires an emissivity en e = 4, a mean free path 
HXuc ~ 1500 kms" 1 , and an intensity nj^e = eHe (HXu e /c) = 0.02. The implied fraction 
of He ii is a^Hen = a rmH e n e / (<7HeC^He^j,He) — 10 _1 A. Therefore, over most of the volume 
in the IGM where A < 1, the Hen fraction is already reduced to less than 10% by the 
low-luminosity sources. Within a gap where the transmitted flux is increased to 20%, we 
need An e — 50 (Fig. 9), increasing the recombination rate by a factor 3 and the mean 
free path by a factor 4 (Figs. 6 and 10), and therefore the intensity by a factor 12, so the 
Hen fraction is reduced to ^g eII — 0.008A. 

Why might the emissivity from quasars be less than is estimated from the observations? 
One possibility is that the He n -ionizing radiation from many quasars is absorbed locally in 
their host halos, or near the nucleus. Only two of the four quasars for which the He n Lyct 
forest has been observed show the strong proximity effect that is expected (Anderson et al. 
1998 and references therein). Among the quasars of lower luminosity that dominate the 
emissivity, a greater fraction of the radiation could be locally absorbed. Another possibility 
is that the value of Qb is higher than predicted by nucleosynthesis for the deuterium 
abundance found by Buries & Tytler (1998). 

Songaila & Cowie (1996) have reported a decrease of the Siiv/Civ ratio with increasing 
redshift for intermediate column density hydrogen Lya absorption lines (see also Boksenberg 
1998). Songaila & Cowie interpreted this evolution as a softening of the spectrum with 
increasing redshift and argued for a rather sudden reionization of Hem at z ~ 3. As 
discussed in §3, reionization is expected to be gradual, and we should therefore expect a 
progressive decrease of the intensity of the Hen -ionizing background with redshift, as the 
Hen photon mean free path decreases. Heap et al. (1999) mentioned that a sudden change 
of the intensity of the He n -ionizing background at z ~ 3 is supported by the sudden 
change in the mean transmitted flux in the Hen Lya spectra. This sudden change in the 
mean transmitted flux may be due to two reasons: first, as seen in Figure 9, even if An e 
decreases smoothly with redshift, the flux decrement reaches unity very rapidly because the 
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voids have a minimum value of their densities A ~ 0.1, as we discussed earlier for the case 
of hydrogen. Second, simple Poisson fluctuations in the number of gaps caused by luminous 
sources may explain the apparently sudden change in the filling factor of these gaps near 
z — 3, since the Hen Lya spectrum has been observed in only four quasars. 

A gradual reionization is also favored by the more recent results of Boksenberg, 
Sargent, & Rauch (1998), where a smooth evolution of the ratios Siiv/Civ and Cn/Crv 
is reported, and also by the findings of Dave et al. (1998) who used the observed Ovi 
absorption to argue that at least 50 percent of the volume of the universe at redshift near 
three is illuminated by a hard spectrum with significant flux well above the Hen Lyman 
edge. 

The reionization of Hen can affect the evolution and spatial distribution of the 
temperature in the IGM (Miralda-Escude & Rees 1994). Each helium atom yields an energy 
input of ~ 54/(7 — l)eV. For 7 = 1.8, this results in a temperature increase of 20000 
K during reionization, which affects the subsequent thermal history of the gas while the 
cooling is long compared to the Hubble time. The equilibrium temperature of photoionized 
gas may also fluctuate substantially owing to the heating by helium photoionization, on 
the scales of the regions of increased Hen ionization around luminous quasars at z ~ 3. 
These temperature fluctuations could affect the evolution of the IGM and the formation of 
low-mass galaxies, because the collapse of the gas is significantly slowed down in halos of 
low velocity dispersion (see e.g., Efstathiou 1992, Thoul & Weinberg 1996). In the likely 
case where the ionizing radiation from a quasar is anisotropic, the present distribution 
of dwarf galaxies might show a large quadrupole moment around the locations of ancient 
quasars, probably occupied today by massive clusters of galaxies. These special galaxy 
fluctuations, which would constitute a new type of non-local galaxy bias unrelated to the 
primordial fluctuations, might be detected in future galaxy redshift surveys on the very 
large scales that Hem regions can reach, where the amplitude of primordial fluctuations is 
greatly decreased. 

6. Discussion and Conclusions 

We have discussed in this paper how the dumpiness of the matter distribution and 
the discreteness of ionizing sources affects the reionization of hydrogen and helium and the 
appearance of the absorption spectra of high-redshift sources. In a clumpy medium, the 
underdense regions filling most of the volume of the universe will be reionized first, while 
the denser regions are gradually ionized later from the outside. As the emissivity increases, 
a balance is maintained with the global recombination rate by increasing the mean free 
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path of ionizing photons, as the size of the neutral, self-shielded regions at high density 
(which are observed as Lyman limit systems) shrinks. 

The first gaps in the hydrogen and helium Gunn-Peterson trough in the absorption 
spectra of high-redshift objects should be caused by the most underdense voids. As the 
redshift is increased, the transmitted flux in Lya spectra decreases rapidly not only because 
the intensity of the ionizing background diminishes if the emissivity from sources is lower, 
but also because the global recombination rate increases, and the voids are less underdense 
as predicted by the gravitational evolution of large-scale structure. We showed that if the 
emissivity does not increase with redshift at z > 4, then the Gunn-Peterson trough should 
be reached at z ~ 6. This redshift is determined by the fact that the most underdense voids 
become optically thick to Lya photons. However, the low-density IGM is likely to have 
been ionized everywhere at this redshift already, and the epoch of overlap of the H n regions 
around individual sources could have occurred at a substantially higher redshift if the 
emissivity is dominated by sources of low luminosity. The rapid increase of the He n Lya 
flux decrement at z ~ 3 can be similarly explained. 

Whether it will be possible to see gaps of transmitted flux in a Lya spectrum across the 
Gunn-Peterson trough, produced by individual H n regions around ionizing sources before 
the major fraction of the volume of the universe is reionized, depends on the luminosity 
of the sources. For a QSO with luminosity 10 46 ergs _1 at z — 6, creating an Hn region of 
radius ~ 70 comoving Mpc, the mean transmitted flux through the Hn region would only 
be 10% (from Figures 5 and 6), and this transmitted flux decreases rapidly for sources of 
lower luminosity. If luminous quasars exist but are short-lived, there should be a first epoch 
of overlap of fossil H n regions (where the size of the H n regions would be smaller by a 
factor (tqeif) 1 / 3 relative to the case of long-lived sources), and a second epoch of overlap of 
active regions, having a size depending only on the luminosity of the quasars. 

Smaller H n regions can on rare occasions produce gaps of transmitted flux when the 
intensity of radiation is much higher than average, when an underdense void that is close 
to the source compared to the size of the Hn region is intersected by the line of sight. 
This should of course happen in the case of the "proximity effect". For an Hn region with 
proper radius smaller than 1 Mpc surrounded by neutral IGM, any residual flux that might 
be transmitted in these circumstances would still be absorbed in the damping wings of the 
absorption by the neutral IGM (see M98). 

Reionization might be complete at a substantially earlier time than the redshift were 
the Gunn-Peterson trough is first encountered, if it is produced by low-luminosity sources. 
We have argued that this is likely to be the case for the Hen reionization, where only a 
modest emissivity from low-luminosity sources is required to reionize the Hen in small 
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Heni regions well before the much larger Hem regions produced by luminous QSO's can 
overlap. The presence of this low-luminosity sources is supported by the presence of a small 
transmitted flux across the Hen Lya spectra at z ~ 3.1 detected by Heap et al. (1999). 

Since reionization occurs outside-in and the dense gas can stay neutral until well after 
the overlap of H n regions, the dumpiness of the IGM does not necessarily increase the mean 
number of times that a baryon will recombine during the reionization epoch; moreover, this 
mean number of recombinations does not just depend on the properties of the IGM, but 
also on the luminosity of the ionizing sources. Thus, for the density distribution suggested 
by numerical simulations and low-luminosity sources, the number of photons per hydrogen 
atom required for the reionization is of order unity. In fact, dumpiness can even reduce the 
required number of photons when a large fraction of baryons are in high density regions 
that do not need to be ionized in order for the mean free path of the photons to increase to 
values much larger than the mean separation between sources. Only at very high redshift, 
when R u ^> 1, recombinations imply a significant increase in the required emissivity for 
reionization. This modest requirement can be met by both faint QSO's and star-forming 
galaxies. The remaining uncertainties to determine which sources reionized the IGM are 
the space density of faint high-redshift QSO's and star-forming galaxies, and the escape 
fraction of Lyman continuum photons from high-redshift galaxies. 

We would like to thank Tom Abel, Piero Madau, and Tom Theuns for helpful 
discussions. We also thank R. Cen and J. P. Ostriker for giving us permission to show in 
Figure 1 the density distribution obtained from their numerical simulation. JM thanks the 
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